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This  paper  presents  the  application  of  Mixed-Integer  Programming  (MIP)  approach  for  solving  the 
security-constrained  daily  hydrothermal  generation  Scheduling  which  takes  into  account  the  inter- 
mittency  and  volatility  of  wind  power  generation,  which  is  called  Security-Constrained  Wind  Hydro- 
thermal  Coordination  (WHTC).  In  restructured  power  systems,  Independent  System  Operators  (ISOs) 
execute  the  Security-Constrained  Unit  Commitment  (SCUC)  program  to  plan  a  secure  and  economical 
hourly  generation  schedule  for  the  daily/weekly-ahead  market.  The  objective  of  security-constrained 
daily  hydrothermal  generation  scheduling  is  to  determine  an  optimum  schedule  of  generating  units  for 
minimizing  the  cost  of  supplying  energy  and  ancillary  services  with  considering  network  security 
constraints.  The  problem  formulation  includes  dynamic  ramp-rate  constraints  for  generation  schedules 
and  reserve  activation,  and  minimum  up-time  and  down-time  of  conventional  units.  Of  particular 
interest  in  this  study  are  considering  more  practical  constraints  and  rigorous  modeling  of  thermal  and 
hydro  units  such  as  prohibited  operating  zones  and  valve  loading  effects.  Furthermore,  for  the  hydro 
plants,  multi  performance  curve  with  spillage  and  time  delay  between  reservoirs  are  considered.  To 
assess  the  efficiency  and  powerful  performance  of  mentioned  method,  a  typical  case  study  based  on 
modified  IEEE-118  bus  system  is  investigated  and  the  results  are  compared  to  each  other  in  different  test 
system. 
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1.  Introduction 

The  electric  energy  securance  issue  has  been  getting  noticed  in  the 
governments  in  whole  of  the  world.  Especially  in  recent  years  that  the 
energy  consuming  has  been  increased  due  to  industry  deploy  and 
dramatic  use  of  electric  devices  instead  of  mechanic  types.  Recently, 
many  of  governments  are  decided  to  invest  in  the  new  resources 
particularly  sustainable  energy  resources  to  providing  the  required 
energy.  The  main  reasons  of  this  plan  are  summarized  as  (a)  as  to 
being  competed  the  fossil  energy  resources  and  raising  its  price  in 
international  energy  market  and  (b)  the  legislation  of  environmental 
laws  such  as  stated  in  the  “Low  Carbon  Transition  Plan",  July  2009,  the 
governments  are  obliged  to  reduce  the  Carbon  emission  to  80%  by 
2050  [1],  Videlicet,  wind  turbine  generators  (WTG)  will  supply  20%  of 
the  US  power  demand  by  2030  [2]  and  European  Union  planned  for 
20%  power  generation  from  renewable  energy  sources  by  2020  [3], 
But,  the  increment  in  penetration  of  wind  power  plants  on  power 
systems  poses  well-documented  technical  challenges  for  the  existing 
power  systems  because  of  uncertainty  production  characteristic. 
Especially,  in  modern  systems  so-called  Smart  Grids,  it  is  necessary 
to  consider  the  unit  commitment  (UC)  and  economic  dispatch 
incorporating  wind  farms  generation. 

Therefore,  uncertainty  of  renewable  energies  must  be  included 
as  an  aim  function  in  system  scheming  daily  and  weekly  means  an 
economic  and  secure  schematization  which  consists  of  the  hydro 
and  thermal  units  incorporating  wind  power  plants.  Daily  Hydro- 
thermal  Generation  Scheduling  (DHGS)  is  a  unit  commitment  and 
economic  dispatch  among  hydro  and  thermal  units  during  a 
scheduling  period  of  time  for  supplying  the  demand  in  minimum 
total  cost.  Also  to  obtain  the  realistic  results,  the  practical  con¬ 
straints  related  to  thermal  plants,  hydroelectric  system  and  elec¬ 
trical  power  system  (i.e.  satisfying  power  system  demand)  should 
be  investigated.  On  the  other  hand,  ISO  coordinates  market 
activities  for  satisfying  the  hourly  load  demand,  limited  fuel  and 
other  resources,  environmental  constraints,  and  transmission 
security  requirements  [4-6], 

Several  literatures  are  studied  on  methods  to  estimation  of 
accurate  wind  speed  and  consequent  power  output  of  wind  farms. 
These  investigations  have  been  based  on  such  foundations  as  fuzzy 
logic  [7],  neural  networks  [8],  and  time  series  [9],  In  this  paper,  the 
time  series  model  of  Auto  Regressive  Integrated  Moving  Average 
(ARIMA)  has  been  used  for  modeling  of  wind  power  generation 
uncertainty  [10], 

Providing  the  Security-Constrained  Hydrothermal  Coordination 
(SCHTC)  with  the  lowest  cost  is  the  nonlinear  problem  with  multiple 
minima.  There  are  some  techniques  such  as  Lagrangian  Relaxation  (LR) 
[4],  Dynamic  Programming  (DP)  [5],  Mixed  Integer  Programming 
(MIP)  [11],  Benders  Decomposition  (BD)  [12]  and  various  intelligent 
techniques  [13-15]  to  find  the  best  solution  for  a  goal  function  with 
complex  and  nonlinear  characteristics  and  heavy  equality  and  inequal¬ 
ity  constraints.  In  this  paper,  the  MIP  approach  is  applied  for  solving 
the  SCHTC  with  volatile  wind  power  generation  problem  which  many 
of  the  practical  constraints  for  detailed  modeling  of  the  all  of  type 
units  are  considered.  Ref.  [16]  develops  a  two-stage  stochastic  mixed- 


quadratic  programming  problem  has  been  proposed  a  new  modeling 
of  the  optimal  unit  commitment  and  sale  bid  to  the  day-ahead  market, 
and  the  optimal  economic  dispatch  of  the  bilateral  contracts  for  all  the 
thermal  and  combined  cycle  units.  This  paper  has  been  used  from  real 
data  of  a  Spanish  generation  company  and  market. 

In  electricity  market,  the  ISO  has  an  important  role  in  manage¬ 
ment  of  the  whole  network  by  possession  of  the  commit  and 
dispatch  of  all  energy  resources  based  to  required  load  with 
meeting  of  security  constraints  and  minimum  cost  [17-19],  Actu¬ 
ally,  the  ISO  determine  an  optimum  schedule  of  generation  units 
with  the  standard  market  design  (SMD)  according  to  security- 
constrained  unit  commitment  (SCUC),  for  minimizing  the  cost 

[20] .  This  is  done  while  in  several  studies  the  GENCO  effort  to 
scheduling  the  generation  of  units  to  obtain  the  maximum  profit 

[21] ,  In  this  paper  the  unit  commitment  and  scheduling  of  system 
units  (hydro,  thermal  and  wind)  is  done  and  analyzed  from  ISO's 
viewpoint  for  having  the  minimum  cost. 

Finally,  in  this  paper  the  both  thermal  and  hydro  subsystems 
are  considered  in  planning  of  ISO  against  of  [20-23]  that  studied 
them  separately  without  considering  of  security  constraints.  Also 
limitations  of  the  ramp  rate,  frequent  switching,  frequent  cycling 
of  hydro  units  in  daily  operations  according  to  mechanics  pro¬ 
blems  and  costs  related  to  start-up  and  wear  and  tear  of  those,  and 
prohibited  operating  zones  (POZs)  in  thermal  units  are  considered 
in  UC  and  ED  [24,25],  In  addition  to  prior  mentions,  the  hydro  and 
thermal  units  operating  characteristics  such  as  valve  point  loading 
and  minimum  up  and  down  time  constraints  are  added  to  the 
proposed  model  for  more  precision  although  may  be  caused  to  a 
higher  nonlinear,  non-smooth  and  non-convex  function  [24].  Also 
to  conducting  the  simulation,  the  load  is  not  constant  and  has 
been  considered  as  a  variable  parameter.  The  load  uncertainty  is 
modeled  based  on  the  error  of  load  prediction  and  needed 
scenarios  for  solving  the  problem  are  generated  by  the  Monte 
Carlo  Simulation  (MCS)  method  and  roulette  wheel  mechanism 
[26,27], 

The  remaining  sections  of  this  paper  are  organized  as  follows. 
In  Section  2,  the  stochastic  model  of  hybrid  wind-hydrothermal 
scheduling  is  formulated  as  a  0/1  mixed-integer  linear  program¬ 
ming  problem.  The  modeling  approaches  used  in  this  paper  to 
represent  the  load  and  wind  power  generation  uncertainties  are 
described  in  Section  3.  Obtained  results  from  case  studies  with 
detailed  discussion  are  presented  in  Section  4,  and  at  the  end  the 
conclusions  are  drawn  in  Section  5. 


2.  MIP  formulation  of  security-constrained  wind 
hydrothermal  coordination 

Mathematically  security-constrained  wind  hydrothermal  coor¬ 
dination  is  a  decision  problem  with  an  objective  to  be  minimized 
with  respect  to  a  series  of  prevailing  equality  and  inequality 
constraints.  There  are  many  parameters,  equality  and  inequality 
constrains  for  reality  optimum  scheduling  of  units  in  power 
system.  These  limitations  must  be  considered  in  the  goal  function 
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Nomenclature 

Indices 

I  thermal  unit  index 

J  hydro  unit  index 

wf  wind  farm  index 

t  time  interval  (hour)  index 

s  scenario  index 

u  uncertain  parameter  index 

Constants 


rj  conversion  factor  equal  to  3.6  x  10  3(Hm3  s/m3  h) 

(1  Hm3  =  106m3) 

0  number  of  periods  of  the  planning  horizon  (24  h) 

Qjj,  t)  minimum  water  discharge  of  unit  j  at  hour  t  (m3/s) 
0(j,  t)  maximum  water  discharge  of  unit  j  at  hour  t  (m3/s) 

Tjj  time  delay  between  reservoir  of  plant  i  and  reservoir 

of  plant  j  (h) 

Aj  shut  down  cost  of  unit  i  ($) 

Aj  start-up  cost  of  unit  j  ($) 

bn(i )  slope  of  block  n  of  fuel  cost  curve  of  unit  i($/MWh) 

bn(j )  slope  of  the  volume  block  n  of  the  reservoir  associated 

to  unit  j  (m3/s/Hm3) 

b„(j)  slope  of  the  block  n  of  the  performance  curve  of  k  unit 
j  (MW/m3/s) 

be„(i )  slope  of  segment  n  in  emission  curve  of  unit  i 
DT(i )  minimum  down  time  of  unit  i  (h) 
et  valve  loading  coefficient 

emin(!)  minimum  emission  generation  of  unit  i  (lbs) 

ECR  emission  group  consists  of  S02,  NOx 

£(p“  j(i))  emission  generation  of  (n  — 1)  th  upper  limit  in  emis¬ 
sion  curve  of  unit  i 

EMmaxO)  maximum  emission  of  unit  i  (lbs) 

EMmix(gr)  maximum  emission  of  group  (lbs) 
ft  valve  loading  coefficient 

FGR  fuel  group  consists  of  coal,  gas  and  oil 

Fmin(i)  minimum  fuel  consumption  of  unit  i  (Mbtu) 

Fmax(i)  maximum  fuel  consumption  of  unit  i  (Mbtu) 

Fmin(gr)  minimum  fuel  consumption  of  group  (Mbtu) 

Fmax(gr)  maximum  fuel  consumption  of  group  (Mbtu) 

F(p“ _,(!))  cost  of  generation  of  (n-l)th  upper  limit  in  fuel  cost 
of  unit  i 

F(j,t,s)  forecasted  natural  water  inflow  of  the  reservoir  asso¬ 
ciated  to  plant  j  in  period  t  (Hm3/h) 

I<A(i )  Cost  of  the  Ath  discrete  interval  of  the  start-up  cost  of 

unit  i  ($/h) 

f°(i)  initial  status  of  unit  i  (0/1) 

L  number  of  variable  head 

M  number  of  prohibited  operation  zones 

NS  number  of  scenario  after  scenario  reduction 
p(s)  probability  of  scenario  s 

p  norm(s)  normalized  probability  of  selected  scenario  s 
JVfSR(i)  maximum  sustained  ramp  rate  (MW/min) 

MU  maximum  number  of  units  that  can  be  on  at  the 
same  time 

NB  number  of  bilateral  contract 

NFS  number  of  discrete  intervals  of  the  start-up  emission 
function 

NE  number  of  discrete  intervals  of  emission  function 
NF  number  of  blocks  of  the  piecewise  linearization  of  the 
variable  cost  function 


pbm(t)  power  capacity  of  bilateral  contract  m  at  hour  t  (MW) 
Pmm(i)  minimum  power  output  of  unit  i  (MW) 

Pmax(Q  maximum  power  output  of  unit  i  (MW) 

Pn(j)  minimum  power  output  of  plant  j  for  performance 
curve  n(MW) 

p(j)  capacity  of  plant  j  (MW) 

pb(i)  lower  limit  of  nth  prohibited  zone  of  unit  i  (MW) 

upper  limit  of  (n  - 1 )  th  prohibited  zone  of  unit  i  (MW) 
Q(j)  minimum  water  discharge  of  hydro  plant  j  if  is  on  (m3/ 

_  s) 

Qn0)  maximum  water  discharge  of  block  n  of  plant  j  (Hm3) 
QSC(i)  quick  start  capacity  of  unit  i  (MW) 

RDLn(i)  ramp  down  limit  for  block  n  (MW) 

RULn(i )  ramp  up  limit  for  block  n  (MW) 
s°(i)  time  periods  of  unit  i  has  been  shut-down  at  the 
beginning  of  the  planning  horizon  (h) 
s(j)  maximum  spillage  of  unit  j  (m3) 

Smax(Q  maximum  hour  unit  i  can  be  off  (h) 

SUE(i)  emission  generated  by  unit  i  when  started  up  (lbs) 

SDE(i)  emission  generated  by  unit  i  when  shut  down  (lbs) 

SDF(i )  fuel  consumption  of  unit  i  when  shut  down  (Mbtu) 
SD(i )  shut-down  ramp  rate  limit  of  unit  i  (MW/h) 

SU(i )  start-up  ramp  rate  limit  of  unit  i  (MW/h) 

UT(i)  minimum  up  time  of  unit  i  (h) 

U°(i)  time  periods  of  unit  i  have  been  on-line  at  the 
beginning  of  the  planning  horizon  (h) 
v0(j)  minimum  content  of  the  reservoir  associated  to  plant  j 

(Hm3) 

v°0)  reservoir  content  at  the  beginning  of  the  study  time 

(Hm3) 

v0(j)  reservoir  content  at  the  end  of  the  study  time  (Hm3) 
v„0)  maximum  content  of  the  reservoir  j  associated  to  nth 
variable  head  (Hm3) 


Variables 

Sn(i,  t)  generation  of  block  n  of  fuel  cost  curve  of  unit  i  at 
hour  t 

/(i,t)  dummy  variable  (h) 

i//  n(i,  t)  generation  of  block  n  of  unit  i  at  hour  t  of  valve  point 
loadings  curve 

y/n(j,  t)  volume  block  n  for  the  reservoir  of  hydro  plant  j  at 
hour  t  (MW) 

B(i,t )  start-up  cost  of  unit  i  at  hour  t  ($) 
bn(i)  slope  of  power  block  n  of  fuel  cost  curve  of  unit  i 
($/MWh) 

bj,0)  slope  of  the  block  n  of  the  performance  curve  1  of 
hydro  plant  j  (MW/m3/s) 

c(i,t )  valve  point  loadings  cost  of  unit  i  at  hour  t  ($) 

F(i,t)  fuel  cost  of  unit  i  at  hour  t  ($) 

Nd(i,t )  non-spinning  reserve  of  a  unit  i  at  hour  t  for  bid  on 
spot  market  when  unit  is  off  (MW) 

Nd(j,t )  non-spinning  reserve  of  a  unit  j  at  hour  t  for  bid  on 
spot  market  when  unit  is  off  (MW) 

Nu(i,t)  non-spinning  reserve  of  a  unit  i  at  hour  t  for  bid  on 
spot  market  when  unit  is  on  (MW) 

Nu(j,t)  non-spinning  reserve  of  a  unit  j  at  hour  t  for  bid  on 
spot  market  when  unit  is  on  (MW) 
p(i,t)  real  power  generation  of  unit  i  at  hour  t  (MW) 

Pmin(h  F,  s),  Pmax(h  F,  s)  lower  and  upper  limit  of  real  power 
generation  of  unit  i  at  hour  t  (MW) 
p(J,t )  real  power  generation  of  unit  j  at  hour  t  (MW) 
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p(j,  t,  s),  p(j,  t,  s )  lower  and  upper  limit  of  real  power  generation 
of  unit  j  at  hour  t  (MW) 

ps(t)  power  for  bid  on  spot  market  at  hour  t  (MW) 

Q0, t)  water  discharge  of  unit  j  at  hour  t  (m3/s) 
qn0,t)  water  discharge  of  block  n  of  unit  j  at  hour  f(m3/s) 
R(i,t)  spinning  reserve  of  a  unit  i  at  hour  t  for  bid  on  spot 
market  (MW) 

R(j,t)  spinning  reserve  of  a  unit  j  at  hour  t  for  bid  on  spot 
market  (MW) 

RDL(p(i,t ))  ramping  down  limit  of  unit  i  at  hour  t  (MW) 
RUL(p(i,t))  ramping  up  limit  of  unit  i  at  hour  t  (MW) 
s(i,t)  time  periods  that  unit  i  has  been  shut-down  at  hour 
t(h) 

s(J,t )  spillage  of  the  reservoir  associated  to  unit  j  at  hour  t 

(m3/s) 

v(J,t)  water  content  of  the  reservoir  associated  to  plant  unit 
j  at  hour  t  (Hm3) 

Binary  variables 

f)n(i,  t,  s)  1  if  block  n  of  fuel  cost  curve  of  unit  i  at  hour  t  selected 
t,s)  1  if  variable  head  n+1  of  unit  j  at  hour  t  selected 
Xn(i,  t,  s )  1  if  power  output  of  unit  i  at  hour  t  has  exceeded  block 
n 

hn(j,t,  s)  1  if  the  water  discharge  of  unit  j  at  hour  t  has  exceeded 
block  n 

/( i, t, s)  1  if  unit  i  is  on-line  at  hour  t 

I(j,t,s)  1  if  hydro  plant  j  is  on-line  at  hour  t 


Id(i,t, s)  1  if  unit  i  at  hour  f  provide  non-spinning  reserve  when 

unit  is  off. 

Idn(i,t, s)  1  if  block  n  of  ramping  down  limit  curve  of  unit  i  at 
hour  t  selected 

Iun(i,t, s)  1  if  block  n  of  ramping  up  limit  curve  of  unit  i  at  hour  f 
selected 

wA(i,  t,s )  1  if  unit  i  is  started-up  at  the  beginning  of  hour  f  and  it 
has  been  offline  for  A  hours 

y(i,t,s)  1  if  unit  i  is  started-up  at  the  beginning  of  hour  t 

y(j,t, s)  1  if  unit  j  is  started-up  at  the  beginning  of  hour  f 

z(i,t,s )  1  if  unit  i  is  shut-down  at  the  beginning  of  hour  t 

z(J,t, s)  1  if  unit  j  is  shut-down  at  the  beginning  of  hour  t 

Sets 

C  set  of  indices  of  the  group  units 

/  set  of  thermal  units 

J  set  of  hydro  units 

WF  wind  farms 

N  set  of  indices  of  the  blocks  of  the  piecewise  lineariza¬ 

tion  of  the  unit  performance  curve. 

T  set  of  indices  of  the  periods  of  the  market  time 

horizon 

S  scenario  numbers 

SP  stochastic  parameters  (w/ indicates  wind  farms) 
i,p,q,r  sets  of  AR1MA  family  model 

A  set  of  the  discrete  intervals  of  the  start-up  cost 

function  for  thermal  units 
Qj  set  of  upstream  reservoirs  of  plant  j. 

N  time  number  of  wind  power  generation  time  series. 


as  security  constraints  included  the  transmission  flow  and  bus 
voltage  constraints  or  the  start-up  limitations  caused  extra  invest¬ 
ment.  In  this  section,  the  objective  function  and  its  different  parts 
of  simulation  is  explained  clearly. 

2.1.  Objective  function 

In  this  work,  the  suggested  objective  function  to  determine  the 
optimal  planning  of  coordination  of  wind,  hydro  and  thermal  units 
over  scenarios  caused  minimum  cost  is  described  as  (1).  It  is 
necessary  to  remind  that  the  load  in  this  paper  is  varied  con¬ 
stantly. 

Costrocal  =  J)  pnorm(s).  cos  t(s) 

S 

cos  t(s)=  £  \  £[F(i,t,s)  +  AZ(i,t,s)  +  B(i,t,s)  +  C(i,t,s)]+  'ZAjy(j,t,s)  l 
'eT0eJ  juj  ) 

(1) 

The  last  term  of  variable  cost  function  of  each  scenario 
represents  the  corresponding  cost  of  hydro  units  regarding  their 
start-up  cost  over  the  given  period  [28],  While  the  first  term 
shows  thermal  operating  cost  including  fuel,  shutdown,  startup 
costs  and  valve  point  loadings  cost.  The  operating  costs  of  wind 
units  are  assumed  to  be  zero.  The  list  of  symbols  is  presented  in 
the  Nomenclature  section. 

When  minimizing  the  total  cost  in  power  systems,  the  total 
generation  of  hydro,  wind  and  thermal  plants  should  be  equal  to 
the  total  system  demand  plus  the  transmission  network  loss.  Due 
to  simplicity,  the  network  loss  is  not  considered  in  this  paper.  This 
equation  is  given  by 

£p(i,t,S)+ Xp(/,t,S)  +  X  p(wf,t,s)  =  PD,t  VteT.VseS  (2) 

is  I  js]  w/eWF 


2.2.  Cascaded-hydro  units'  model 

According  to  previous  sections,  the  hydro  power  plants  have 
the  start-up  limitation.  Because  of  the  unnecessary  commitments, 
loss  of  water  during  start-up  period,  wear  and  tear  of  the  windings 
and  mechanical  equipment,  the  start-up  limitation  must  be 
considered  in  planning.  In  this  regard,  the  relation  among  gener¬ 
ated  power,  water  discharge,  and  multi  performance  curves  of 
hydro  power  plant  is  shown  in  Fig.  1.  The  proposed  model  takes 
into  account  head  dependent  reservoirs  with  MIP  formulations  as 
well  as  hydro  plants  connected  both  in  parallel  and  series  shown 
in  Fig.  2.  The  number  of  the  heads  for  power  plants  is  considered 
to  be  L  in  the  following  formulations  represented  in  Fig.  1. 

Performance  curves  based  on  the  water  volume  can  be  given  by 
(3)— (6)  show  that  the  volume  of  each  hydro  power  plant  must  be 
bigger  than  the  minimum  content  of  that  hydro  plant.  Eqs.  (4)  and 
(5)  determine  the  right  head  proportional  to  volume.  Combination 


Generation  (MW  ) 


Fig.  1.  Three-dimensional  piecewise  linear  non-concave  unit  performance  curve 
for  hydro  plant  j  at  hour  t. 
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Fig.  2.  Hydraulic  topology  of  the  river  basin. 


of  0-1  for  these  binary  variables  is  prevented  by  Eq.  (6).  /?„  (J,t,s) 
equals  to  1  when  (n+l)th  head  is  used  or  water  volume  of  the 
reservoir  is  greater  than  Vn(J).  In  other  words,  these  equations 
choose  the  right  curve  for  head  according  to  the  content  level. 

v(j,t,S)>v0(j)  Vj  eJ,Vt  eT,Vs  eS  (3) 

v(j,  t,S)<vL(j)/3L_1(j,  t,s) 

+  Z  vn_iO')[^n_20,  t,S)]  Vj  ej,  VteT, VseS  (4) 

n  =  2 

v(j,  (,S)>Vn(j)/),_,(j,  t,s ) 

+  Z  Vn_20’)[/in_20',t)-y3n-i0',f,5)]  VjeJ,  VteT,  VseS  (5) 

n  =  3 

/?i(j,  t,s)>/32(j.  t,s)>  ...  >f)L_i(j,t,S)  VjeJ,  VteT,  VseS  (6) 


2.2.3.  Piecewise  linear  formulation  of  the  variable  head  multi 
performance  cwves 

In  most  papers  the  effect  of  the  variation  of  the  head  has  been 
neglected  and  head  is  assumed  fix.  The  proposed  M1LP  model  take 
into  consider  head  dependent  reservoirs  with  a  piecewise  linear 
approximation.  After  determination  of  each  unit  in  this  subsection, 
the  linear  relationship  between  the  generated  powers,  the  dis¬ 
charged  water  and  variable  head  for  corresponding  performance 
curve  can  be  presented  as  follows: 

p(j,  t,s)-pk(j)l(j,  t,  s)-  £  qn(j,t,s)bkn(j ) 

neN 


- plj )  (k— 1)—  z /?„(/,  t,S)+  Z/?n0',t,s) 
n  =  1  n  =  k 

VjeJ,VteT,VseS,1<k<L 


<0 


(7) 


p(j,LS)-pk(j)l(j,t,S)-  2  <7n0>  t,S)bn(j) 
neN 


(k— 1)—  Z  /J„(j,t,S)  +  Z  /^n0’ t,  s) 

n = 1  n=k 


+Pd ) 

VjeJ, VteT, VseS,1<k<L 


>0 


(8) 


where  p(j,t,s)  is  the  generated  power  by  the  hydro  unit  j  at  hour  t, 
Pk(j)  is  the  minimum  power  generation  of  the  head  k  which  is 
determined  by  pn(j,t,s).  Also,  P(j)  is  the  capacity  of  hydro  unit  j,  and 
<JnC/.t,s)  is  the  water  discharge  of  the  block  n.  Finally,  bnk(j )  is  the 
slope  of  the  block  n  of  the  variable  head  k  of  hydro  unit  j. 


2.2.2.  Water  discharge  limits 

The  water  discharge  limit  and  initial  and  final  volume  and 
water  balance  equations  same  as  the  equations  presented  in  [22]; 
however,  in  the  proposed  WHTC  model  the  spillage  effect  is  also 
considered  which  is  inspired  by  [29], 

2.3.  Thermal  units'  model 

In  this  subsection  deal  with  the  linearization  of  all  nonlinear 
equations  related  to  thermal  units  to  represent  the  linear  model 
for  thermal  units. 

2.3.3.  Linear  fuel  cost  function  considering  POZ 

The  thermal  units  have  some  prohibited  operating  zone  due  to 
steam  valve  operation  or  vibration  in  its  shaft  bearing  and  some 
faults  in  the  machines  or  their  accessories  such  as  pumps  or 
boilers,  etc.  It  has  been  shown  that  the  input-output  curve  of 
thermal  units  has  some  valve  points.  These  valve  points  generate 
some  prohibited  zones,  so  the  input-output  curve  of  thermal  unit 
will  be  discontinuous.  The  quadratic  production  cost  function  can 
be  accurately  approximated  by  a  set  of  piecewise  blocks.  So,  the 
linear  formulation  of  power  generation  for  the  ith  thermal  unit 


considering  M  POZs  is  modeled  as  follows: 

M+l 

F(i,  t,  s)  =  £  [fin(i,t,s)F(pun_^i )) 

n  =  1 

+ bn(i)Sn(i,  t,s)]  Vie/,  Wt  eT,Ws  g  S  (9) 

M+l 

P(i,t,s)=  Z  [Pn-l(i)/?„(/,t,S) 

n  =  1 

+<5„(i,f,s)]  Vi  el,  VteT,  VseS  (10) 

Sn(i,t,s)>  0  n  =  l,2,...,M+  L,VieI,  VteT,  VseS  (11) 

s„(i,  t,  S)  <  [Pn(t)— (i)]/3n(t,  t,  s)  n=  1,2,  ...,M+  1, 

Vi  el,  VteT,VseS  (12) 

M+l 

2  t, s )  =  /(i,  t,s)  Vi  g  I,Vt  g  T,Vs  g  S  (13) 

n  =  1 

PniUUs)  G  {0,1}  n=  1,2,  +  1,  Vie/,  VteT,  VseS  (14) 

where,  p“(i)  =  pmin(i)  and  p^+1  (i)  =  pmax(i). 

2.3.2.  Cost  of  valve  point  loadings 


In  practical,  the  generators  have  the  multiple  valves  in  steam 
turbines.  The  opening  and  closing  of  these  valves  caused  to 
maintain  the  active  power  balance.  However  it  adds  the  ripples 
in  the  cost  function  [30],  To  calculate  and  consider  the  valve-point 
effects,  sinusoidal  functions  are  added  to  the  quadratic  cost 
functions  as  follows  [31-33]: 

Fi(Pi)  =  ai  +  biPi+ciP2  +  abs(eisin(fi(Pimjn-Pi)))  Vi  el  (15) 

where  e,  and  f  are  the  coefficients  of  generator  ith  considering 
valve  point  loading  effect. 

But  the  added  sinus  term  to  cost  function  cause  to  having  a 
nonlinear  and  non-convex  function  that  can  be  obstacle  in  MIP 
model.  In  this  paper  a  linear  formulation  is  proposed  to  overcome 
this  problem  as  shown  in  Fig.  3. 

C(i,t,s)  =  '^{  +2  Z  Wau+i (i,t,s)-i//4n+4(i,t,s)] 

K  [  n  =  0 

+(2-v/2)  2  [i//4n+2(i,  t,  S)-i//4n+3(i,  t,  s)]  1  Vi  el,  VteT,  VseS 

n  =  0 


(16) 
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Fig.  3.  Linear  approximation  absolute  sinus  function  valve  point  loading. 


where,  y/n  is  power  generated  by  nth  block  in  sth  scenario. 
p(i,  t,  S)  =  pmin(i)I(i,  t,  S)  +  j  [v/'4n+l(i>t>s)  +  V 4n+2&t’S) 

n  =  0  L 

+V/4n+3(,>t=S)  +  V/4n+4(i>t’S)]  Vi  S  I,  V  t  6  T,  VS  6  S  (17) 

-^-X\(i,t,S)<y/^(i,t,S)<-^-I(i,t,S)  Vie/,  Vf  eT,  VseS  (18) 

^vT„(i,  t, s)  <  f, s)  <  t.s)  Vi  e  I,  V t  e  T,  Vs  e  S,  n  =  2, 3, . . xf 

(19) 

Xn(i,t,s)e  {0,1}  VieI,VteT,VseS,n  =  i,2,...,Xi  (20) 

where, 

/<,•  =floorlfi(pmax(i)-pmin(i)/yr)]  and  xt  =floor[4fi(pmax(i)-pmin(i)/jr)] 

Constraint  (17)  states  that  power  output  of  unit  i  at  hour  t  is  the 
sum  of  minimum  power  output  if  unit  is  on,  plus  the  power 
generated  in  each  block.  Constraint  (18)  limits  the  power  gener¬ 
ated  in  first  block,  this  power  should  be  greater  than  zero  and 
smaller  than  or  equal  to  n/4fi  that  is  “power  length”  of  each  block. 
In  this  constraint,  I(i,t,s)  is  used  ensuring  that  power  generated  of 
first  block  equal  to  zero,  if  unit  i  is  off  at  hour  t.  To  limit  the 
generated  power  in  each  block  xn(h  t,s)  is  introduced  in  con¬ 
straints  (18)-(20).  This  binary  variable  equal  one  if  the  generated 
power  of  unit  i  at  hour  t  has  exceeded  block  n. 

2.3.3.  Dynamic  ramping  up/down  limit 

Variations  of  load  have  been  responded  by  generators  by 
decreasing  and  increasing  of  generating  power.  But  the  rate  of 
the  power  output  changing  is  limited  by  two  constraints.  The 
increasing  and  decreasing  constraint  is  defined  by  up  and  down 
ramp  rate  that  are  specified  for  each  generator.  In  this  work  a 
dynamic  ramp  rate  is  used,  as  the  ramp  rate  is  a  function  of  the 
generating  power.  Also  in  this  method  the  M  POZ  are  considered. 
The  formulation  related  to  dynamic  ramp  rate  is  given  as 

M+l 

RUL(p(i,t,s))  =  2  RULn(i)fin{i,t,s)  Vi  e  /,V  t  eT,Vs  e  S  (21) 

n  =  1 
M+l 

RDL(p(i,t,s))  =  Z  RDL„(i)/3n(i,t,s )  VieI,VteT,VseS  (22) 

n  =  1 

Constraints  (21)  and  (22)  indicate  dynamic  ramp  limit  with  M 
POZs.  The  binary  variables  are  used  to  indicate  which  operating 
zone  has  been  selected.  Another  formulation  for  dynamic  ramp 
rate  is  proposed  in  [34], 

2.3.4.  Stair-wise  formulation  of  the  time-dependent  startup  cost 
function 

Starting  up  generating  unit  costs  money  in  addition  to  the 
running  cost  considered  in  economic  dispatch.  This  cost  is  mod¬ 
eled  as  a  nonlinear  (exponential)  function  of  the  number  of  hours 
a  unit  has  been  off.  In  [35,36]  two  different  methods  have  been 


suggested  for  time  dependent  linearization  of  this  function.  In  this 
paper  is  used  of  formulation  proposed  in  [35], 

2.3.5.  Minimum  up-time  (MUT)  and  minimum  down-time  (MDT) 
Another  major  source  of  complication  in  system  planning  is  the 
nonlinear  minimum  up/down  time  constraints  which  has  to  be 
checked.  In  this  study  a  linear  formulation  for  these  constraints  is 
represented  as  follows: 

B(i) 

Z  [l-/(Us)]  =  0  vie/  (23) 

t=  l 

r, 

Z  I(i,T,s)>  UT(i)y(i,t,s)  T-,  =  t  +  UT(i)-l, 

T  =  t 

Vt  =  B(i)  +  1, ....  @—UT(i)  +  1  (24) 

0 

z  [/(i,T,s)-y(i,t,s)]  >0  VieI,Vt  =  0—UTi  +  2,...,&  (25) 

T  =  t 

C(i) 

Z  l(i,t,S)=  0  Vie/  (26) 

t  =  l 

z  [1— /(i,T,s)]  >DT(i)Zi(t,s)  T2  =  t  +  DT(i)- 1. 

T  =  t 

Vt  =  C(i)  +  l,...,<9-DT(i)+l  (27) 

0 

z  [1-l(i,T,s)-z(i,t,s)]>0  Vi  el,  Vt  =  0-DT(i)  +  2, ....  0  (28) 

r  —  t 

Once  the  operation  of  a  unit  begins,  the  operational  status  should 
be  maintained  as  least  for  a  certain  time  (MUT).  Eq.  (23)  deter¬ 
mines  the  initial  status  of  units  so  unit  will  be  on  if  B(i)  is  less  than 
UT  (i).  The  MUT  limit  for  consecutive  periods  and  for  the  last  hours 
is  satisfied  by  (24)  and  (25)  respectively.  Similarly,  if  the  operation 
of  a  unit  is  terminated,  the  un-operational  status  should  be 
maintained  at  least  for  a  certain  time  (MDT)  [36],  Constraints 
(26)-(28)  enforce  the  MDT  limit. 

The  detailed  formulation  of  other  hydro  units  constraints  such 
as  initial  and  final  volume  [22],  water  balance  [29]  and  operating 
services  [37]  considered  in  the  problem,  are  represented  in  the 
mentioned  references. 

2.4.  Security  constraints 

The  security  constraints  required  in  security-constrained  wind 
hydrothermal  coordination  problem,  are  obtained  based  on  DC 
power  flow  model  due  to  it  has  more  accuracy  than  the  linear 
power  flow  model.  In  this  model  the  physical  flow  in  a  grid  is 
controlled  by  Kirchoffs  current  law  (KCL)  and  Kirchoffs  voltage 
law  (KVL),  while  in  linear  power  flow  model  only  the  KCL  is 
considered.  In  DC  power  flow  model,  the  transmission  constraints 
are  considered  as  linear  constraints.  So  the  practical  constraints 
related  to  thermal  plants,  hydroelectric  system  are  expressed 
by  piecewise  linear  approximation,  therefore  the  security-con- 
strained  WHTC  is  a  Mixed  Integer  Linear  Programming  (MILP) 
problem.  The  related  formulation  to  the  DC  power  flow  is  given  as 


Fit  =  ^r0is~^ir)  Vb,VteT  (30) 

Transmission  flow  limits  in  the  base  case: 

-FrX  ^  f  it  ^  f  l?aX  Vl.VteT  (31) 
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3.  Wind  power  generation  and  demand  forecast  uncertainty 

One  of  the  most  important  issues  considered  recently  by 
researchers  of  restructure  power  systems  widely,  is  the  uncertainty 
in  operation  of  power  systems.  Uncertainties  such  as  uncertainty 
due  to  system  load  prediction  error  or  uncertainty  due  to  unex¬ 
pected  departure  of  network  equipment  (i.e.  production  units  and 
lines).  One  another  these  uncertainties  related  to  the  error  in 
prediction  of  the  accurate  value  of  wind  farms  output,  especially 
which  penetration  of  this  energy  resource  type  in  power  grids  are 
increasing. 

3.1.  Wind  power  uncertainty  modeling 

The  wind  power  is  a  non-predictable  complete  accurately.  This 
future  is  rooted  from  the  wind  nature.  Therefore  the  primary 
problem  for  incorporating  wind  power  plants  in  power  system  is 
the  uncertainty  of  them  that  decrease  the  reliability  of  grid. 
Several  investigations  have  looked  at  the  prediction  of  wind  speed 
for  use  in  determining  the  available  wind  power  [7,8], 

In  literatures  that  study  on  the  wind  farms  uncertainty,  the 
wind  power  plant  is  modeled  on  two  method,  wind  speed-based 
approaches  [38-42],  and  wind  power-based  approach  [10],  In  the 
wind  speed-based  approaches,  the  wind  speed  must  be  measured 
accurately  and  the  generating  power  be  calculated  with  model  of 
wind-turbine  that  usually  is  unavailable.  While  in  power-based 
approach  does  not  require  the  model  of  wind  farm  due  to  the 
output  of  them  is  measured  and  recorded  directly,  so  the  whole  of 
wind  power  data  with  its  details  is  available.  Also  in  this  method 
the  error  of  prediction  is  smaller  due  to  the  amount  of  generated 
power  by  a  wind-turbine  is  proportional  by  the  cube  of  wind 
speed,  so  an  small  error  (e)  in  the  wind  prediction  produce  a  more 
big  error  (e3)  in  expected  power  output  than  when  the  wind 
power  is  measured  directly. 

Although  there  are  other  methods  for  consideration  of  volatility 
of  the  wind  power  generation  like  as  Monte  Carlo  simulation,  and 
Latin  hypercube  sampling  method  [43],  in  this  paper,  inspired  by 
[44],  an  Auto-Regressive  Integrated  Moving  Average  (ARIMA) 
process  is  used  to  model  wind  power  generation  uncertainty 
directly  as  a  wind  power-based  stochastic  approach.  In  the 
following,  the  ARIMA  process  is  described  clearly. 

3.1.1.  ARIMA  model  of  wind  power 

An  ARIMA  ( p,d,q )  model  of  the  non-stationary  random  process 
W(f)  is  depicted  as  [10]: 

l-.Z  <PiD^(a-DYW(t)  =  d0+  (V  £  6>,D!'jn(t)  (32) 

where  [qtf]  is  the  Auto-Regressive  (AR)  coefficients;  {#,}  is  the 
moving  average  (MA)  coefficients;  n(t)  is  a  white  Gaussian  process 
with  zero  mean  and  variance  the  parameter  0o  refers  to  the 
deterministic  trend  term  when  r  >  0. 

The  coefficients  of  the  ARIMA  model  can  be  estimated  by 
different  approaches,  e.g.  the  Yule-Walker  estimator,  the  least 
square  estimator,  and  the  maximum  likelihood  estimator  [45],  For 
specification  of  the  ARIMA's  family  there  are  two  model  selection 
criteria  in  various  stochastic  analyses  that  the  Akaike's  Information 
Criterion  (AIC)  [45]  and  Bayesian  Information  Criterion  (BIC)  is 
more  common  than  other  types.  The  ARIMA(0,1,1)  model  that  is 
appropriate  for  the  wind  power  time  series  is  described  as  Eq.  (33). 

(l-D)Wo(t)  =  (l-d^D)a(t) 

WP(t)  =  Wl(t)  for  t=  1, ....  N 


In  this  paper,  the  ARIMA(0,1,1)  are  used  in  wind-hydrothermal 
coordination  for  generating  of  different  wind  power  generation 
scenarios  and  modeling  of  random  trend  of  the  wind  farm. 

The  relation  between  the  input  wind  power  and  the  output 
electric  power  relies  on  several  factors,  such  as  the  efficiencies  of 
generator,  wind  rotor,  gearbox,  and  inverter,  depending  on  what 
kind  of  turbine  under  investigated.  For  a  generic  wind  turbine, 
some  researchers  [46,47]  used  a  simplified  model  to  characterize 
the  relation  between  the  wind  power  generation  (WPG)  and  wind 
speed  [48]: 

w  =  0  for  v<Vj  and  v>v0  (34) 

v— vr 

w  =  wr -  for  v,<v<v0  (35) 

Vr-Vj 

w  =  wr  for  vr  <  v  <  v0  (36) 

Eqs.  (34)-(36)  show  the  power  output  of  wind  power  plants  is 
limited  by  wind  speed,  while  the  ARIMA  model  cannot  consider 
these  limitations.  Thus  these  margins  has  been  considered  for 
output  wind  farms  and  added  to  the  standard  ARIMA  model  to 
include  the  upper  and  lower  bounds  of  the  WPG  as  follows: 

fWmax  W0(t)>W 

max 

W0(t)  =  <  Wo(t)  Wmin<W0(t)<W  max  (37) 

[  Wmin  Wo(t)  <  Wmin 

where  Wmax  and  Wmin  denote  the  upper  and  lower  bounds  of  the 
square  root  of  the  wind  power  output,  respectively.  The  block 
diagram  of  the  limited  ARIMA(0,1,1 )  is  shown  in  Fig.  4. 

3.2.  Load  uncertainty  modeling 

In  this  paper,  the  load  uncertainty  is  considered  as  the  error  of 
load  prediction.  Thus,  the  probability  distribution  function  of  load 
prediction  error  must  be  derived  based  on  the  records  of  earlier 
loads.  It  should  be  mentioned;  here  it  is  assumed  that  the  load 
forecasting  process  as  one  of  the  tasks  of  the  ISO,  before  the 
implementation  market  has  been  done.  Another  noticeable  point 
is  that,  considering  each  of  the  loads  connected  to  the  network  as 
an  independent  variable  will  be  caused  to  increasing  of  the 
number  of  random  variables  and  the  problem  will  be  more 
complex.  To  solve  this  problem,  the  network  load  is  considered 
as  a  random  variable.  It  is  natural  that  for  getting  the  value  of  each 
of  the  loads  connected  to  the  network,  it  is  sufficient  to  be  divided 
based  on  the  ratio  of  each  bus  load  to  the  total  load  in  the  base 
state.  Fig.  5  is  an  example  of  the  probability  distribution  function 
of  the  prediction  error  of  network  load  and  with  its  discretization. 

As  is  seen  in  Fig.  5,  seven  different  load  intervals  centered  mean 
error  of  zero  (base  state)  is  considered,  so  that  distance  of  between 
the  different  levels  of  load  forecast  error  equal  to  the  standard 
deviation  (SD)  of  load  prediction  error. 

The  important  issue  is  how  to  model  the  random  load  level  as  a 
random  variable.  To  implement  this  is  used  from  the  roulette 
wheel  mechanism  [49,50].  For  this  aim,  first  is  applied  the 
possibility  of  different  levels  of  load  prediction  error  based  on 
the  unit  load  so  that  the  total  probability  is  equal  to  one.  Then, 
the  distance  between  range  [0,1]  is  covered  based  on  the  normal 


(33) 


Fig.  4.  Block  diagrams  of  ARIMA(0,l,l)-type  wind  power  time  series  models  [21], 
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Fig.  5.  Typical  discretization  of  the  probability  distribution  of  the  load  error 

[26,27]. 
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Fig.  6.  Roulette  wheel  mechanism  for  the  normalized  probabilities  of  the  load 
forecast  levels. 


probability  of  forecast  error  levels  of  load.  Naturally,  whatever  the 
probability  of  prediction  error  levels  of  load  is  more;  it  will  occupy 
more  space  of  the  roulette  wheel.  To  clarify  this  mechanism,  the 
roulette  wheel  in  Fig.  6  is  illustrated. 

After  forming  of  the  roulette  wheel,  a  random  number  in  the 
range  [0,1]  is  generated.  The  number  produced  will  be  in  specified 
ranges  on  the  roulette  wheel  are  associated  with  different  levels  of 
load  forecast  error,  this  means  that  the  load  level  related  to  this 
interval  are  selected  for  the  desired  scenario.  Thus,  for  generating 
of  each  new  scenario,  the  roulette  wheel  mechanism  will  be  used 
to  model  the  behavior  of  random  load. 

A  higher  number  of  scenarios  results  in  a  better  modeling  of 
uncertainties  but  with  the  cost  of  higher  computation  burden. 
Accordingly,  after  the  production  scenarios,  it  may  be  a  very  low 
probability  of  some  scenarios  and  will  only  cause  to  increasing  of 
time  and  volume  calculations.  Therefore,  it  is  necessary  to  elim¬ 
inate  the  Low-value  scenarios  from  the  set  of  produced  scenarios 
and  decrease  the  number  of  scenarios.  In  this  paper,  the  strategy 
adopted  for  this  purpose  is  the  removing  scenarios  with  less  than 
probability  and  the  similar  scenarios  [51].  This  reduction  process 
results  a  model  of  uncertainties  with  more  less  accuracy  and 
negligible  error  with  the  cost  of  higher  computation  burden. 

3.3.  Scenario  aggregation 

The  idea  of  stochastic  security-constrained  WHTC  is  to  con¬ 
struct  or  sample  possible  options  for  uncertain  circumstances, 
solve  the  deterministic  security-constrained  wind  hydrothermal 
coordination  problem  for  the  possible  options,  and  select  a  good 
combination  of  the  outcomes  to  represent  the  stochastic  solution. 
Two  methods  are  usually  considered  for  scenario  aggregation  (the 
combination  of  the  scenario  outcomes)  of  the  stochastic  security- 
constrained  wind  hydrothermal  coordination  [52]:  (1)  using  mini¬ 
mum  value  of  the  random  variables  over  the  scenarios  and 
(2)  applying  weighted-average  (expected  value)  of  variables  over 
the  all  scenarios.  Using  minimum  value  of  daily  operating  cost  may 


lead  to  pseudo-scheduling  with  optimistic  operating  costs  in  the 
planning  horizon.  Because  it  may  be,  minimum  daily  operating 
cost  be  related  to  the  scenario  with  minimum  occurrence  prob¬ 
ability  and  has  a  significant  difference  ratio  to  other  scenarios 
determined  daily  operating  cost  of  problem.  So,  can  be  inferred 
that  using  of  first  method  is  a  very  conservative  method  and  not 
economical.  While  in  second  method  videlicet  weighted-average 
or  expected  value,  scenarios  with  more  probability  will  have  a 
more  contribution  on  specification  of  final  value  of  daily  operating 
cost.  It  seems  which  second  method  is  more  reasonable  of 
economic  view  point.  For  this  reason,  in  this  paper  the  second 
method  is  used  for  aggregation  of  different  scenarios  result  paper 
to  determine  total  daily  operating  cost  of  problem. 

So,  in  the  proposed  stochastic  optimization  framework,  the 
later  method  (expected  value)  is  considered.  In  this  way,  the 
solutions  obtained  from  different  scenarios  are  aggregated  based 
on  the  probability  laws  to  yield  a  single  solution,  describing  the 
most  probable  outcome  of  the  power  system  based  on  the 
evaluated  scenarios,  considered  as  the  result  of  the  proposed 
stochastic  security-constrained  wind  hydrothermal  coordination 
framework. 

As  stated,  the  MCS  method  is  implemented  to  simulate  random 
characteristics  of  power  systems  load  and  then  the  scenario 
aggregation  technology  is  used  to  convert  the  stochastic  variables 
of  the  stochastic  security-constrained  wind  hydrothermal  coordi¬ 
nation  problem  into  deterministic  ones.  A  major  advantage  of 
scenario  aggregation  technique  is  that  not  only  individual  scenar¬ 
ios  become  simple  to  interpret  but  also  the  underlying  problem 
structure  is  preserved.  After  running  the  proposed  security- 
constrained  wind  hydrothermal  coordination  scheme  for  the 
accepted  scenarios  resulted  from  the  scenario  reduction,  the 
results  are  aggregated  according  to  the  probability  of  scenarios 
to  get  the  expected  results  of  the  formulation  of  hybrid  wind- 
hydrothermal  scheduling  considering  uncertainties. 

The  aggregation  is  done  for  the  scenario  dependent  decision 
variables  /(i,t,s),  I(j,t, s),  F(i,t, s),  p(i,t,s),  p(j,t,s),R(i,t,s),  R(j,t,s )  of  the 
optimization  problem.  The  aggregation  is  done  as 

NS 

f  =  I  Prm-fs  Vs  eS  (38) 

S  =  1 

where,  /  is  the  variable  that  is  aggregated  and  fs  is  the  variable 
value  at  scenario  s.  It  is  noted  that  the  objective  function  of  the 
proposed  formulation  for  security-constrained  WHTC  problem  in 
Eq.  (1)  is  also  an  aggregation  of  the  objective  function  values  of  the 
scenarios. 

In  the  present  paper,  two  main  categories  of  power  system 
uncertainties  (load  and  wind  power  uncertainties)  are  considered 
in  security-constrained  wind  hydrothermal  coordination  program 
that  execute  by  ISO  to  plan  a  secure  and  economical  hybrid  wind- 
hydrothermal  generation  scheduling. 

When  the  number  of  variation  and  constraint  increased  by 
increasing  the  size  of  system,  the  solution  time  will  increased  and 
sometimes  solver  cannot  solve  the  problem,  but  there  are  several 
method  to  overcome  this  obstacle.  Benders  decomposition  approach 
[53]  is  useful  method  for  this  situation,  also  parallel  processing  may  be 
used  for  solving  this  problem,  for  solving  method  with  non  linear 
constraint,  linearization  maybe  handsome. 


4.  Case  studies 

A  modified  IEEE  118-bus  test  system  is  used  to  test  the  proposed 
algorithm  for  the  security  constrained  day-ahead  wind-hydrothermal 
power  scheduling  [54].  In  addition  to  hydro  and  thermal  units,  two 
wind  farms  are  assumed  which  are  supervised  financially  under  ISO 
operating  scheme.  Random  behavior  of  wind  power  generation  is 
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Table  1  Table  2 

Results  of  the  scenarios  of  the  stochastic  security-constrained  WHTC  framework.  Scenario  10:  WHTC  with  security. 


No. 

Scenario  no. 

Normalized 

probability 

Total  wind 
power  (MW) 

Daily  operating  cost  ($) 

Daily  operating  cost:731830.98  $ 
Hours  (0-24) 

1 

1 

0.788 

1572.175 

755704 

1-3 

1000000000000000000000000 

2 

26 

0.030 

1238.43 

839811 

4 

1111111111111111111111111 

3 

28 

0.013 

1639.317 

802221 

5 

1111111111111111111111110 

4 

39 

0.010 

1557.42 

742289 

6 

1000000000000000000000000 

5 

40 

0.008 

1393.175 

861473 

7 

1111111111111111100000000 

6 

45 

0.009 

1466.166 

761167 

8-10 

1000000000000000000000000 

7 

52 

0.010 

1257.439 

865732 

11 

1111111111111111111111111 

8 

56 

0.015 

1497.414 

867398 

12-18 

1000000000000000000000000 

9 

78 

0.005 

1305.838 

801530 

19 

1111111111111111111111110 

10 

79 

0.033 

1551.793 

731831 

20-21 

1111111111111111111111111 

11 

82 

0.013 

1250.692 

854589 

22-23 

1000000000000000000000000 

12 

87 

0.003 

1565.795 

770844 

24-25 

1111111111111111111111111 

13 

97 

0.006 

1445.347 

775077 

26 

1000000000000000000000000 

14 

101 

0.009 

1353.77 

850058 

27-29 

1111111111111111111111111 

15 

105 

0.024 

1504.765 

777797 

30-38 

1000000000000000000000000 

16 

133 

0.006 

1504.793 

789590 

39 

1111111111111100000000000 

17 

145 

0.005 

1408.149 

804971 

40 

1100000000000000000000000 

18 

158 

0.001 

1640.804 

826181 

41-43 

1000000000000000000000000 

19 

164 

0.006 

1542.076 

765741 

44-45 

1111111111111111111111111 

20 

194 

0.006 

1232.111 

895407 

46-51 

1000000000000000000000000 

modeled  by  ARIMA( 0,1,1 )  time  series  model  as  described  in  Section  3.1. 
Two  wind  farms  are  separately  modeled  to  simulate  correlation 
concerns  of  the  stochastic  wind  power  behavior  in  different  areas  of 
the  network.  Auto  and  cross  correlation  coefficients  related  to  two 
applied  wind  farms  are  considered  as  Appendix  I.  As  mentioned  in 
Section  3,  the  proposed  MIP  model  of  security-constrained  wind 
hydrothermal  coordination  includes  both  uncertainties  of  wind  power 
and  load.  The  practical  constraints  of  thermal  and  hydro  generation 
units  that  detailed  their  modeling  are  mentioned  in  Section  2  and 
additional  system-wide  constraints  such  as  fuel  constraints  and 


53-54  1000000000000000000000000 

Hydro 1  1111111111111111111111110 

Hydro2  1111111111111111100000000 

Hydro3  1111111111111111111111110 

Hydro4  1111111111111111111111110 

Hydro5  1111111111111111111111110 

Hydro6  1111000000000001111111100 

Hydro7  1111111111111111111111110 

Hydro8  1111111111111111111111111 


Table  3 


emission  limits  [36,38,55]  and  spinning  and  operating  reserve  require¬ 
ments  [56]  are  considered  in  the  stochastic  optimization  framework. 

The  proposed  model  is  implemented  on  a  Pentium  IV,  3  GHz 
with  1GB  RAM  using  MILP  solver  CPLEX  9.0  in  the  GAMS  environ¬ 
ment  [57]. 

With  the  use  of  ARIMA  method  for  modeling  the  stochastic 
behavior  of  wind  power  generation  and  roulette  wheel  mechanism 
to  model  uncertainty  of  load,  200  scenarios,  including  daily  wind 
power  time  series  and  daily  load  profiles,  are  generated.  It  imposes  a 
high  computational  burden  to  solve  the  security-constrained  WHTC 
problem  for  all  of  these  scenarios.  So,  the  set  of  generated  scenarios 
(200  wind  power  generation  scenarios/  daily  load  profile)  is  reduced 
using  the  scenario  reduction  technique.  The  generated  similar  scenar¬ 
ios  and  scenarios  with  probability  lower  than  0.003  are  discarded. 
Number  of  remaining  scenarios  after  scenario  reduction  is  equal  to  20, 
which  results  in  200/20=10  filtering  ratio.  So,  the  scenario  reduction 
technique  significantly  reduces  the  computation  burden  of  the  pro¬ 
posed  stochastic  security  constrained  WHTC  framework.  At  the  same 
time,  the  most  probable  and  dissimilar  scenarios  are  retained  while 
maintaining  a  good  approximation  of  the  uncertain  behavior  of  these 
uncertainty  resources.  For  the  remaining  set  of  scenarios,  the  proposed 
stochastic  MILP  model  of  security-constrained  wind  hydrothermal 
coordination  is  run.  Selected  scenarios,  their  normalized  probability 
and  total  daily  power  generation  of  wind  farms  (generated  power  by 
WF1  and  WF2)  are  presented  in  Table  1. 

Case  1.  Security  constrained  wind-hydrothermal  coordination 

In  state  that  security  constraints  are  considered,  the  minimum 
and  maximum  daily  operating  cost  are  related  to  scenario  10  and 
scenario  20  with  731831$  and  895407.38  $  respectively.  The 
commitment  schedules  for  these  scenarios  are  shown  in 


Scenario  20:  WHTC  with  security. 


Daily  operating  cost:895407.39  $ 
Hours  (0-24) 

1 

1000000000000000000000000 

2-3 

1000000000000000000011111 

4-5 

1111111111111111111111111 

6 

1000000000000000000000000 

7 

1000000000000000000001111 

8-10 

1000000000000000000000000 

11 

1111111111111111111111111 

12-18 

1000000000000000000000000 

19-21 

1111111111111111111111111 

22-23 

1000000000000000000000000 

24-25 

1111111111111111111111111 

26 

1000000000000000000000000 

27-29 

1111111111111111111111111 

30-35 

1000000000000000000000000 

36 

1000000000000000000001111 

37-33 

1000000000000000000000000 

39 

1111000000000000000000000 

40 

1111111111111111111111111 

41-42 

1000000000000000000000000 

43 

1100000000000000000000000 

44 

1111000000000000000000000 

45 

1111111111111111111111111 

46-51 

1000000000000000000000000 

52 

1111111111111111111111111 

53-54 

1000000000000000000000000 

Hydro  1 

1111001111111111111111110 

Hydro2 

1111111111111111111111110 

Hydro3 

1111001111111111111111110 

Hydro4 

1111111111111111111111010 

Hydro5 

1111111111111111111111111 

Hydro6 

1000000010000000011111110 

Hydro7 

1000001101110111100101011 

Hydro8 

1111111111111111111111111 
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Tables  2  and  3,  in  which  1  and  0  represent  ON/OFF  states  of  units 
at  different  hours,  and  hour  0  represents  the  initial  condition. 

According  to  simulation  of  scenario  10  and  20,  the  following 
results  can  be  obtained: 

In  scenario  10,  the  thermal  units  2  and  3  which  have  the  high 
production  cost  is  remained  off  on  all  period  planning.  While  in 
scenario  20,  both  these  units  in  the  final  of  the  day  hours  (from 
hour  20  to  24)  have  been  forced  to  turn  on  and  connected  to 
the  grid. 

In  scenario  10,  the  thermal  plant  40  is  not  committed  at  hour 
2-24,  while  in  scenario  20  this  plant  is  committed.  In  scenario  10, 
the  thermal  plant  44  is  ON  at  total  of  24  h  while  in  scenario  20  it  is 
ON  only  at  hour  1  -3.  The  expensive  units  of  plant  52  is  committed 
at  hour  1-13  in  scenario  10  while  it  will  be  ON  at  total  of  24  h  in 
scenario  20.  Also  the  expensive  plant  7  and  39  has  been  ON  more 
hours  in  scenario  10  than  scenario  20.  So  expensive  plant  7  is 
connected  to  grid  from  hour  1  to  16  while  in  scenario  20  it  has 
been  committed  only  at  hour  21-24. 

The  hydro  units  6  and  7  in  scenario  10  are  committed  more 
hours  than  scenario  20. 

The  whole  generating  power  of  hydro  and  thermal  plants  in  the 
scenarios  10  and  20  are  15107.94/54218.871  MW  and  14656.002/ 
62380.987  MW  respectively.  As  can  be  seen,  in  scenario  with 
lower  production  cost,  i.e.  scenario  10  the  amount  of  hydro  plants 
production  is  more.  Also  the  ratio  of  the  hydro  generating  power 
to  total  of  power  generated  with  hydro  and  thermal  plants  is 
21.79%  and  19.02%  in  scenario  10  and  20  respectively. 

In  scenario  10,  the  thermal  plants  generate  the  54218.871  MW 
power  in  all  period  of  time  and  the  average  of  generating  power 
with  those  is  2259.119  MW.  While  the  maximum  generating 
power  of  thermal  power  is  3369.614  MW  in  hour  15.  It  is 
noteworthy  that  variation  of  thermal  plants  generating  power  is 
less  compared  with  variation  of  hydro  plants  generating  power 
given  in  next  paragraph  ((3369.614-2259.119)/2259.119=49.16%). 
The  reason  of  these  low  variations  is  the  limitation  of  thermal 
plants  power  variations  and  due  to  this  constraint,  the  thermal 
units  do  not  follow  the  rapid  load  variations. 

In  scenario  10,  the  hydro  units  in  all  period  of  time  produce  the 
15107.94  MW  that  its  average  is  the  629.497  MW  at  each  hour. 
The  minimum  generating  power  at  hour  14  is  104.415  MW  and 
the  maximum  of  hydro  units  generating  power  at  hour  15  is 
1393.809  MW.  The  variations  of  hydro  units  generating  power  are 
very  high  and  pursue  the  load  variations  in  instantaneous  market. 
Also  the  amount  of  whole  S02  and  NOX  gases  produced  during 
24  h  has  been  80697.95  and  41878.83  in  scenario  10  and  101391.11 
and  50755.54  in  scenario  20  respectively.  This  indicates  that 
despite  production  cost  in  scenario  10  is  less  than  scenario  20, 
the  amount  of  pollution  gases  generation  is  less  too.  Also  it  should 
be  noted  that  the  minimum  pollution  gases  (S02  and  NOX) 
generation  related  to  scenario  10. 

Case  2.  Wind-hydrothermal  coordination  without  considering 
of  security  constraints 

Again  each  20  scenarios  have  been  simulated  without  con¬ 
sidering  security  constraints  and  their  power  production  cost  has 
been  compared  with  each  other.  It  can  be  seen  that  with  not 
considering  of  the  security  constraints  the  minimum  and  max¬ 
imum  of  production  cost  daily  related  to  the  scenario  1  and  20 
with  698609.36  $  and  811174.28  $  respectively.  The  commitment 
schedule  of  scenario  1  is  shown  in  Table  4.  Also  it  is  obvious  which 
considering  of  security  constraints  will  be  caused  to  increasing  the 
power  production  cost  daily  compared  to  when  they  are  neglected. 
The  main  reason  is  that  if  prior  due  to  satisfying  security 
constraints  some  cheap  units  cannot  be  ON,  now  in  this  state 
can  be  ON.  With  comparing  of  scenario  1  (scenario  with  lowest 


Table  4 

Scenario  1 :  WHTC  without  security. 


Daily  operating  cost:  698636.98$ 
Hours  (0-24) 

1-4 

1000000000000000000000000 

5 

1111111111111111111111111 

6-10 

1000000000000000000000000 

11 

1111111111111111111111111 

12-19 

1000000000000000000000000 

20-21 

1111111111111111111111111 

22-23 

1000000000000000000000000 

24-25 

1111111111111111111111111 

26 

1000000000000000000000000 

27-29 

1111111111111111111111111 

30-38 

1000000000000000000000000 

39-40 

1111111111111111111111111 

41-44 

1000000000000000000000000 

45 

1111111111111111111111111 

46-54 

1000000000000000000000000 

Hydro  1 

1111111111111111111111110 

Hydro2 

1111111111111111111111110 

Hydro3 

1111111111111111111111110 

Hydro4 

1111111111111111111111110 

Hydro5 

1111111111111111111111111 

Hydro6 

1000000000000000111110000 

Hydro7 

1111111111111111111110000 

Hydro8 

1111111111111111111111111 

Table  5 

Congested/uncongested  status  of  violated  lines  in 
scenario  1  (with  minimum  cost)  of  Case  2. 

Hours  (0-24) 

L37 

0000000000000000011010000 

L41 

0000000000000000000110000 

L124 

0000000010001001011100000 

L126 

0000000000001001000100000 

L127 

0000000000001000000000000 

L134 

0000000010001001011110000 

L136 

0000000010001001011110000 

LI  37 

0000000000001001000100000 

L138 

0000000000001000000100000 

L139 

0000000000001000000100000 

L142 

0000000000000000111110000 

cost  related  to  state  without  considering  of  the  security  con¬ 
straints)  and  scenario  10  (scenario  with  lowest  cost  related  to 
state  with  considering  of  the  security  constraints),  it  is  seen  that 
the  both  inexpensive  units  4  and  5  are  committed  in  scenario  10, 
while  in  scenario  1  only  unit  5  maintains  its  ON  state.  Also  it  is 
obtained  that  the  both  expensive  thermal  units  7  and  52  in 
scenario  10  to  supplying  load  and  satisfying  of  security  constraints 
has  been  forced  to  be  ON  at  hour  1-16  and  hour  1-13  respectively, 
whereas  in  scenario  1  these  unit  all  24  h  has  not  been  committed. 
Compared  scenario  1  with  scenario  10,  thermal  unit  40  that  was 
OFF  in  all  of  24  h  has  been  ON  and  unit  39  is  committed  at  the 
more  hours  to  compensate  for  the  reduced  supply  for  decommit¬ 
ting  of  thermal  units  19,  44  and  52  and  to  satisfy  physical 
constraints  in  scenario  1. 

With  comparing  of  scenario  20  (scenario  with  highest  cost) 
with  and  without  considering  of  the  security  constraints,  it  is  seen 
that  the  economical  unit  4  and  both  inexpensive  units  19  and  52 
had  been  forced  to  be  ON  over  a  day  to  supplying  load  and 
satisfying  of  security  constraints,  have  been  OFF.  Instead  of  these 
decommitted  units,  inexpensive  unit  10  that  was  OFF  before  will 
be  ON  at  total  of  24  h.  Also  expensive  thermal  units  2,  3  and  7  and 
units  34  and  44  that  were  ON  only  in  the  first  or  last  hours  of  the 


736 


M.  Karami  et  al.  /  Renewable  and  Sustainable  Energy  Reviews  28  (2013)  726-737 


Table  6 

Congested/uncongested  status  of  violated  lines  in 
scenario  20  (with  maximum  cost)  of  Case  2. 


Hours  (0-24) 


L30  0000000000000000000111000 
L37  0000000100000000000101010 
L41  1000001111000000001101111 
L54  0000000000000000000001000 
LI  24  0000001101000000000001011 
LI  26  0000000001000000000000000 
LI  27  0000000000000000000000000 
LI 34  0000001101000000000101010 
L136  0000001101000000000101010 
LI 37  0000000000000000000000000 
L138  0000000000000000000000000 
L139  0000000000000000000000000 
L142  0000000000000000001101110 


Table  7 

Scenarios  with  minimum/maximum  daily  operating  cost  ($)  in  3  cases. 


With  security  constraints  Without  security  constraints 


Case  1  Case  2  Case  3  Case  1  Case  2  Case  3 


Scenario  10  731831  -  - 

Scenario  20  895407.38  - 

Scenario  1  -  - 

Scenario  15  -  -  991005.67 


811174.28 

698609.36 


Table  S 

Aggregated  results  of  the  stochastic  security-constrained  WHTC  framework  in 
3  cases. 


Case  1 

Case  2 

Case  3 

Daily  operating  cost  ($) 

788596.71 

730831.36 

1086309.34 

day,  and  also  without  considering  security  constraints  during 
these  hours  they  have  been  OFF  and  decommitted  over  a  day. 
For  compensating  of  reduced  supply  due  to  decommitting  these 
5  units,  unit  39  that  was  ON  only  in  first  3  h  of  day  it  have  been 
forced  to  be  ON  over  a  day  to  supplying  load. 

When  in  this  case,  we  calculate  the  WHTC  solution  by  excluding 
transmission  and  voltage  constraints,  transmission  flow  violations  in 
scenario  1  and  scenario  20  occur  in  different  lines  and  hours. 
Tables  5  and  6  show  all  of  transmission  flow  violations  on  congested 
lines  for  these  scenarios  without  considering  security  constraints,  in 
which  1  and  0  represent  congested/uncongested  status  of  lines  at 
different  hours,  and  hour  0  represents  the  initial  condition. 

Results  show  that  branch  41  is  very  congested.  This  branch  has 
80  MW  of  capacity  and  that  is  not  sufficient  to  transmit  less 
expensive  generation  from  the  right-hand  side  of  the  system  to 
the  left-hand  side  and  the  its’  transmission  flow  in  Case  3  is  near 
or  at  the  capacity  limit  at  hours  1,  2,  8,  9, 10, 11, 12, 14, 15, 16, 18  and 
19.  Also  scenario  20  has  less  congested  line  than  scenario  1,  for 
example  congested  lines  127, 137, 138,139  in  scenario  1  do  not  have 
any  violated  flow  in  scenario  20. 

Case  3.  Evaluation  of  effect  of  hydro  units  on 
the  security-constrained  WHTC  problem 

Once  again,  each  20  scenarios  have  been  simulated  and 
investigated.  But  this  time  with  eliminating  of  hydro  plants  from 
problem  and  then  their  power  production  cost  that  related  to  the 
thermal  and  wind  units,  are  compared  with  each  other.  In  this  case 


the  minimum  and  maximum  production  cost  daily  are  related  to 
the  scenario  15  with  991005.67  $  and  scenario  20  with  1204288.30 
$  respectively.  Because  of  not  produce  power  by  hydro  units  which 
are  inexpensive  units,  it  can  be  seen  that  production  cost  in  each 
scenario  is  higher  than  two  previous  cases  included  hydro  units. 
For  quick  reference,  the  three  case  tests  are  briefly  described  in 
Table  7.  Also,  the  aggregated  results  of  the  stochastic  security- 
constrained  wind  hydrothermal  coordination  framework,  accord¬ 
ing  to  the  scenario  aggregation  procedure,  are  given  in  Table  8. 
Using  the  proposed  stochastic  framework,  all  20  accepted  scenar¬ 
ios  contribute  into  determining  the  security-constrained  wind 
hydrothermal  coordination  results  according  to  their  probability 
values. 


5.  Conclusion 

In  this  paper,  the  thorough  and  comprehensive  M1P  formula¬ 
tions  for  solving  security-constrained  wind  hydrothermal  coordi¬ 
nation  problem  at  day-ahead  are  proposed.  The  objective  is  to 
minimize  the  daily  operating  cost  of  ISO.  The  proposed  framework 
consists  of  many  practical  constraints  such  as  prohibited  operating 
zones,  valve  loading  effects,  crew  constraints,  dynamic  ramp  rate, 
and  fuel  and  emission  limitations  for  thermal  units.  Furthermore, 
for  the  hydro  plants,  multi  performance  curve  with  spillage  and 
time  delay  between  reservoirs  are  considered.  We  also  consider 
transmission  network  constraints.  The  main  feature  of  the  pro¬ 
posed  framework  refers  to  the  linear  nature  of  the  formulations 
which  is  very  important  for  application  of  the  model  in  the  large 
scale  and  the  real  size  power  system.  Therefore,  the  proposed 
scheme  is  practical  to  generate  appropriate  information  for  ISOs. 
Roulette  wheel  mechanism  and  Monte  Carlo  Simulation  (MCS)  are 
employed  for  random  scenario  generation  wherein  the  stochastic 
optimization  problem  is  converted  into  its  respective  deterministic 
equivalents.  Results  of  the  implementing  the  proposed  model 
show  that  the  solution  time  of  the  problem  is  rationale.  The 
research  work  under  way  in  order  to  reduce  the  computational 
burden  due  to  increasing  the  problem  size,  two  approaches  (and 
their  combinations)  can  be  taken.  These  are  the  use  of  decom¬ 
position  approaches  such  as  benders  decomposition  and  parallel 
computing  used  for  each  scenario. 
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Appendix  A 

Auto  and  cross  correlation  coefficients  related  to  two  imple¬ 
mented  wind  farms  [10]: 
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